Spatial mapping of cell types across the mouse brain (2/3) - cell2location

This notebook demonstrates how to use the cell2location model for mapping a single cell reference cell types onto a spatial transcriptomic dataset. Here, we use a 10X single nucleus RNA-sequencing (snRNAseq) and Visium spatial transcriptomic data generated from adjacent tissue sections of the mouse brain (Kleshchevnikov et al., BioRxiv 2020).

The first step of our model (#2 in Fig 1, tutorial 1/3) is to estimate reference cell type signatures from scRNA-seq profiles, for example as obtained using conventional clustering to identify cell types and subpopulations followed by estimation of average cluster gene expression profiles (Suppl. Methods, Section 2, Fig S1). Cell2location implements this estimation step based on Negative Binomial regression, which allows to robustly combine data across technologies and batches (Suppl. Methods, Section 2).

Figure 1.
Figure 1. Overview of the spatial mapping approach and the workflow which are enabled by cell2location. From left to right: Single-cell RNA-seq and spatial transcriptomics profiles are generated from the same tissue (1). Cell2location takes reference cell type signatures derived from scRNA-seq and spatial transcriptomics data as input (2, 3). The model then decomposes spatially resolved multi-cell RNA counts matrices into the reference signatures, thereby establishing a spatial mapping of cell types (4).

In the second step covered by this notebook (#4 in Fig 1), cell2location decomposes mRNA counts in spatial transcriptomic data using these reference signatures, thereby estimating the relative and absolute abundance of each cell type at each spatial location (Suppl. Methods, Section 1, Fig S1).

Outline

The cell2location workflow consists of three sections:

I. Estimating reference expression signatures of cell types (1/3)

II. Spatially mapping cell types (2/3, this notebook):

  1. Loading Visium data, QC
  2. Visualise Visium data in spatial 2D and UMAP coordinates
  3. Load reference cell type signature from snRNA-seq data and show UMAP of cells
  4. Cell2location model description and analysis pipeline, Evaluating training

III. Results and downstream analysis (3/3)

Loading packages and setting up GPU

First, we need to load the relevant packages and tell cell2location to use the GPU. cell2location is written in pymc3 language for probabilistic modelling that uses a deep learning library called theano for heavy computations. While the package works on both GPU and CPU, using the GPU significantly shortens the computation time for 10X Visium datasets. Using the CPU is more feasible for smaller datasets with fewer spatial locations (e.g. Nanostring WTA technology).

Tips on initializing GPU THEANO_FLAGS='force_device=True' forces the package to use GPU. Pay attention to error messages that might indicate theano failed to initalise the GPU. E.g. failure to use cuDNN will lead to significant slowdown.

Above you should see a message similar to this confirming that theano started using the GPU:

/lib/python3.7/site-packages/theano/gpuarray/dnn.py:184: UserWarning: Your cuDNN version is more recent than Theano. If you encounter problems, try updating Theano or downgrading cuDNN to a version >= v5 and <= v7.
  warnings.warn("Your cuDNN version is more recent than "
Using cuDNN version 7605 on context None
Mapped name None to device cuda0: Tesla V100-SXM2-32GB (0000:89:00.0)

Do not forget to change device=cuda0 to your available GPU id. Use device=cuda / device=cuda0 if you have just one locally or if you are requesting one GPU via HPC cluster job. You can see availlable GPU by openning a terminal in jupyter and running nvidia-smi.

1. Loading Visium data

In this tutorial, we use a paired Visium and snRNAseq reference dataset of the mouse brain (i.e. generated from adjacent tissue sections). There are two biological replicates and several tissue sections from each brain, totalling 5 10X visium samples.

First, we need to download and unzip spatial data, as well as download estimated signatures of reference cell types, from our data portal:

Now, let's read the spatial Visium data from the 10X Space Ranger output and examine several QC plots. Here, we load the our Visium mouse brain experiments (i.e. sections) and corresponding histology images into a single anndata object adata.

Note! Mitochondia-encoded genes (gene names start with prefix mt- or MT-) are irrelevant for spatial mapping because their expression represents technical artifacts in the single cell and nucleus data rather than biological abundance of mitochondria. Yet these genes compose 15-40% of mRNA in each location. Hence, to avoid mapping artifacts we strongly recommend removing mitochondrial genes.

Look at QC metrics

Now let's look at QC: total number of counts and total number of genes per location across Visium experiments.

2. Visualise Visium data in spatial 2D and UMAP coordinates

Visualising data in spatial coordinates with scanpy

Next, we show how to plot these QC values over the histology image using standard scanpy tools

Find regions determined by mLTO / LTi

Add counts matrix as adata.raw

Construct and examine UMAP of locations

Now we apply the standard scanpy processing pipeline to the spatial Visium data to show experiment to experiment variability in the data. Importatly, this workflow will show the extent of batch differences in your data (cell2location works on samples jointly, see below).

In this mouse brain dataset, only a few regions should be different between sections because we are using 2 samples from biological replicates sectioned at a slightly different location along the anterior-posterior axis in the mouse brain. We see general alighnment of locations from both experiments and some mismatches, but as you will see in the downstream analysis notebook most of the differences between experiments here come from batch effect, which cell2location can account for.

3. Load reference cell type signature from snRNA-seq data and show UMAP of cells

Next, we load the pre-processed snRNAseq reference anndata object that contains estimated expression signatures of reference cell types (see notebook 1/3).

Export reference expression signatures of cell types:

Quick look at the cell type composition in our reference data in UMAP coordinates (UMAP representation was generated using a standard scanpy workflow, see notebook 1/3).

Run NNLS-autogenes

4. Automatic discrete region identification by clustering cell densities

We find regions by clustering spots based on inferred molecule contributions of each cell type. We use leiden clustering that incorporates both the similarity of spots in cell locations and in their proximity, by including both when computing the KNN graph. Results are saved in adata_vis.obs['leiden']. Since the clustering is done jointly, the cluster identities match between sections.

# Compute average abundance of each region cluster from cell2location.cluster_averages.cluster_averages import get_cluster_averages_df from cell2location.plt.plot_heatmap import clustermap ct_aver = get_cluster_averages_df(X=adata_vis.obs[['q05_spot_factors' + i for i in adata_vis.uns['mod']['fact_names']]], cluster_col=adata_vis.obs["region_cluster"]) ct_aver.index = adata_vis.uns['mod']['fact_names'] #ct_aver.columns = ['region_' + c for c in ct_aver.columns] # normalise to get 10% of each cell type in each location ct_aver = (ct_aver.T / ct_aver.max(1)).T with mpl.rc_context({'axes.facecolor': 'white'}): clustermap(ct_aver, fun_type='dotplot',#.iloc[linkage['row_ord'], :], cluster_rows=True, cluster_cols=True, figure_size=(12, 25)) plt.savefig(f"{scvi_run_name}/region_clusters_vs_cell_type.pdf", bbox_inches='tight', facecolor='white')

6. Identify co-occuring cell type combinations using matrix factorisation

Here we use inferred cell densities as input to non-negative matrix factorisation to identify co-occuring cell type combinations.

results_folder+r['run_name']from cell2location import run_colocation res_dict, sp_data = run_colocation( adata_vis, model_name='CoLocatedGroupsSklearnNMF', verbose=False, return_all=True, train_args={'n_fact': np.arange(8, 35), 'n_iter': 20000, 'sample_name_col': 'sample', 'mode': 'normal', 'n_type': 'restart', 'n_restarts': 5}, model_kwargs={'init': 'random', 'random_state': 0, 'nmf_kwd_args': {'tol': 0.001}}, posterior_args={}, export_args={'path': results_folder+r['run_name']+'/CoLocatedComb_tol_0001/', 'run_name_suffix': ''}) from cell2location import run_colocation res_dict, sp_data = run_colocation( adata_vis, model_name='CoLocatedGroupsSklearnNMF', verbose=False, return_all=True, train_args={'n_fact': np.arange(8, 35), 'n_iter': 20000, 'sample_name_col': 'sample', 'mode': 'normal', 'n_type': 'restart', 'n_restarts': 5}, model_kwargs={'init': 'random', 'random_state': 0, 'nmf_kwd_args': {'tol': 0.0001}}, posterior_args={}, export_args={'path': results_folder+r['run_name']+'/CoLocatedComb/', 'run_name_suffix': ''})

Modules and their versions used for this analysis